Random Packings of Frictionless Particles 
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We study random packings of frictionless particles at T = 0. The packing fraction where the 
pressure becomes nonzero is the same as the jamming threshold, where the static shear modulus 
becomes nonzero. The distribution of threshold packing fractions narrows and its peak approaches 
random close-packing as the system size increases. For packing fractions within the peak, there is 
no self-averaging, leading to exponential decay of the interparticle force distribution. 
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A system jams when it develops a yield stress or ex- 
tremely long stress relaxation time in a disordered state 
101 . Different control parameters can be varied to induce 
jamming, such as the temperature T, the applied shear 
stress a, or the packing fraction (f), as shown in the phase 
diagram inset to Fig. |l](a) Q. Such a phase diagram 
might apply e.g. to supercooled liquids, granular materi- 
als, foams and suspensions. For the diagram to be useful, 
there should be a common physical origin for jamming 
independent of the control parameter varied. Previously, 
it was shown ||^ that a peak in the distribution of inter- 
particle normal forces, P{F), signifies the development 
of a yield stress in a variety of systems Q , implying that 
the jamming phase diagram is a useful concept. 

There is a special point on the jamming phase diagram, 
marked J in Fig. |^(a), for repulsive, finite-range poten- 
tials. This point, at zero temperature and zero shear 
stress, represents the onset of jamming with increasing 
packing fraction. Static granular packings must neces- 
sarily lie near this point because they are effectively at 
T = (the thermal energy is much smaller than the en- 
ergy needed to lift a grain by its own diameter) and the 
particles are hard, so it is difficult to compress packings 
further into the jammed region. Experimentally, P{F) 
for granular packings has a remarkably robust form 
not only does it have a peak, but it also has an exponen- 
tially decreasing tail at large F. Numerous simulations 
find the same form for P{F) ^j^. The persistence of 
the exponential tail, independent of the potential used, 
is surprising. 

In this letter, we examine configurations created by 
quenching systems from high temperature to T = (and 
also cr = 0) near the onset of jamming (point J in the 
jamming diagram). We find that different configurations 
have different properties, even for arbitrarily large system 
sizes, so that self-averaging is not observed. However, 
the range of packing fractions over which self-averaging 
is not observed narrows with increasing system size. We 
will show that one consequence of non-self-averaging is 
an exponential decay of P{F) at large F. 
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FIG. 1. (a) Pressure p vs. — where (j>c is the threshold 
packing fraction for a given state. Circles (squares) corre- 
spond to Of = 2(5/2) in 2D and diamonds (triangles) corre- 
spond to a = 2(5/2) in 3D and iV = 1024(512) in 2D{3D). 
Inset: Jamming phase diagram showing onset of jamming at 
point J. (b) Shear stress G vs. (j) — (jic. The shear strain 
is applied in the a;-direction, the strain gradient is in the 
y-direction, and the xy component of the stress tensor is mea- 
sured, (c) Number of overlaps per particle Z — Zc vs. (f> — 4>c. 



To create configurations near point J we start with ran- 
dom configurations (i.e. T — oo) and, as with inherent 
structures ||^ , quench infinitely rapidly to T = at fixed 
(f> using conjugate gradient (CG) energy minimization M. 
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We study 50 : 50 binary mixtures of N frictionless parti- 
cles with a diameter ratio 1.4 [|l^,|l^ that interact pair- 
wise via repulsive potentials: Vab{r) = (e/a) (l — r/aab)" 
for r < (Tab and Vabir) = for r > (Jab-^ where a,b la- 
bel particles and Uab — i^a + o'fc)/2. We study harmonic 
(a = 2) and Hertzian {a = 5/2) potentials in 2D and 
3D. The total potential energy is F = if no particles 
overlap. Energy is measured in units of e and length in 
units of the small particle diameter ai . 

We classify each final configuration as either over- 
lapped {V ^ 0) or non-overlapped = 0). Overlapping 
configurations must have a nonzero pressure, p, while 
non-overlapping ones have p = 0. Is a configuration that 
has p > necessarily jammed with a zero- frequency shear 
modulus G > 0? To answer this question, we study states 
close to overlap threshold as a function of packing frac- 
tion. We find the overlap threshold, 0c j for each con- 
figuration by compressing a non-overlapped state while 
measuring p. Throughout the compression, after each 
small increment in </> (Jc/) < 10~^), we use CG to bring 
the state to the lowest energy attainable without cross- 
ing any barriers. (To check that no barriers were crossed, 
we reproduced our results using ten times smaller incre- 
ments in 0). This procedure allows us to measure zero- 
frequency properties of the system. 

Different states have different values of the overlap 
threshold, Nevertheless, when we plot pressure -p 
versus (j) — (j)c (Fig- 0(a)), the results for different config- 
urations collapse on a single curve. This holds for both 
harmonic (a — 2) and Hertzian (a — 5/2) potentials in 
2D and 3D. We find p = po{(f> — (f>c)'^ , where po is only 
weakly dependent on dimension and /3 = 1.0 for har- 
monic and /3 = 1.5 for Hertzian potentials, independent 
of dimension. This is consistent with previous results at 
larger — 0c [i|,|l2|,OI ■ states with p > are 

jammed, we calculate the shear modulus, G, by apply- 
ing a small step strain and measuring the infinite-time 
response by minimizing the energy using CG. The re- 
sponse is linear for sufficiently small strains. Fig. |l](b) 
shows that G = Go((/) — 4>c)^ , where 7 = 0.5 for a = 2 
[||, and 7 = 1.0 for a = 5/2, in 2D and 31?. To our res- 
olution, which is better than 10^^, we find that G and p 
vanish at the same packing fraction (/)c- This implies that 
the onset of jamming, as defined in the first paragraph, 
coincides with the onset of overlap [p^ . 

When a configuration jams at 0c, the number of over- 
laps per particle, Z, jumps from zero to a threshold value 
Zc- Above 0CJ Z increases as Z — Zc = Zo{ip — (j)c)'^, as 
shown in Fig. ^(c). For both harmonic and Hertzian po- 
tentials, we find C = 0.5 in both 2D and 3D |l2). (If 
there are no zero-energy modes, Zc = 2d for frictionless 
spheres in d dimensions 0|. Usually w 5% of particles 
are "rattlers" that do not overlap with any neighbors. If 
we remove these, we find Zc = 2d.) 

By studying the onset of jamming in repulsive sys- 
tems at T = 0, we have a criterion for whether a state is 



jammed or not {i.e. F > or y = 0). If there is a pack- 
ing with V — 0, then it is an equilibrium configuration. 
The state that maintains V ^ Q up to the highest is 
therefore the one that, when compressed infinitesimally 
above this point, becomes the zero-temperature equiva- 
lent to the ideal glass. We have found that the properties 
shown in Fig. ^depend only on — 0c- This suggests that 
this behavior is the same as for the ideal glass. 
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FIG. 2. (a) The probability fj of finding a jammed state vs. 
for Of = 2 in 3-D and several system sizes, TV. (b) The prob- 
ability distribution Pj{4'c) of finding a jamming threshold (j>c 
for systems considered in (a), (c) Full-width at half-maximum 
w of Pj[(j>c) VS. A'', (d) The deviation in peak position 0o of 
Pj{(j}c) from its asymptotic value 0* vs. A'^. Data for (c) and 
(d) are for systems considered in Fig. ^ except a — 5/2 in 2D. 

Although Fig. |l| shows scaling behavior above 0c, the 
value of 0c varies from state to state. How much can 0c 
vary? Our protocol, of starting with random configura- 
tions and quenching infinitely rapidly to T = should, in 
principle, allow us to sample all of phase space. For each 
0, we measure the fraction of initial states that lead to 
jammed states at T = 0. Fig. ^(a) shows the probability, 
/j(0), of finding a jammed state for different system sizes 
in 3D for the harmonic potential. For each N , we have 
differentiated /j with respect to to obtain the proba- 
bility distribution -Pj(0c) of finding a jamming threshold 
of 0c (see Fig. |](b)). System size has a large effect. For 
each N, we characterize -Pj(0c) by its full width at half 
maximum, w, and its peak position, 0o. Fig. 0(c) shows 
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that for N > 10, w (X N where uj « 0.55. Similarly, 
Fig. ||(d) shows that 0o approaches its asymptotic value 
(j)* as 0* - (/>o oc iV-^ where 6 w 0.7 for 2D and « 0.5 
for 3D. We find 0* = 0.842 in 2D and 0* = 0.648 in 3L», 



unlikely that a jammed system can exist with pockets 
containing more than a few rattlers. 



close to well-known values of random close-packing |15|. 

For finite N, the peak 0o of the distribution Pj(0c) cor- 
responds to the largest number of initial states that lead 
to final jammed configurations (not the largest number 
of distinct final jammed states). As such, <po is a mea- 
sure of the largest fraction of phase space that leads to 
the onset of jamming for a given N. (f>* thus represents 
where the jamming threshold is maximally random in 
the N —f oo limit. We find the same limiting value of 
(j)* for Hertzian and harmonic potentials, suggesting that 
(j)* is not sensitive to the potential. We can approach 
jammed hard-sphere packings by noting that states up 
to the jamming threshold are accessible to hard spheres. 
If we repeat the measurement of <j)* for potentials with 
progressively harder repulsions, we can approach (but 
not attain) the hard-sphere limit. This suggests a way to 
measure the maximally random jammed packing fraction 
for hard spheres | |T^ ]. 

It is well known that spherical granular materials can 
exist over a 15% spread of packing fractions ranging from 
0.55 to 0.64 The width of our distribution Pj{(f>c) 

for isotropic packings of frictionless soft spheres is a max- 
imum near iV = 10 and becomes arbitrarily small in the 
N ^ oo limit. Only for TV « 10 can we find states that 
are jammed at packing fractions as low as 0.55. Thus, the 
experimental difference between loose- and close-packing 
cannot be explained simply by considering allowed config- 
urations of frictionless spheres; the value of loose-packing 
is not a purely geometrical quantity. 

We note that there are other protocols for finding the 
jamming threshold at T = 0. We have also generated 
configurations by cooling slowly from equilibrium ther- 
mal states to T = 0. In that case, we find that Pj{4>c) 
is shifted to higher values of (f)c, with values of 0o that 
are less than 1% higher than for the first protocol, but 
it is difficult to determine whether the difference persists 
when N —^ oo. 

We have shown that as iV — + cx), the range of pack- 
ing fractions over which systems can jam becomes ar- 
bitrarily narrow. Therefore, one might suppose that the 
fact that different configurations jam at different packing 
fractions might become irrelevant in the large N limit. 
This is not the case because within the range over which 
both jammed and unjammed configurations exist, there 
is no self-averaging. For example, if a configuration is un- 
jammed, with V — 0, every subset of that configuration 
has V — 0, a,s well. Less obvious is that for jammed con- 
figurations, all subsets of more than a few particles will 
likely also contain overlaps. This is found numerically, 
and stems from the fact that on average, each jammed 
particle must have at least 2c? overlapping contacts, each 
of which also has 2d contacts. This constraint makes it 
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FIG. 3. (a) Force distribution P{F/{F)) obtained by scal- 
ing F by the average force {F) of each configuration using the 
harmonic potential in 3-D with A'^ = 1024. The solid line is a 
Gaussian fit to the high-force tail, (b) Same as (a) except F 
is scaled by {{F)) averaged over all configurations. 

The distribution of interparticle normal forces, P{F), 
illustrates the absence of self-averaging. Simulations of 
static packings 1^,0 have shown that at large F, P{F) 
falls off exponentially even for harmonic and Hertzian 
potentials. An argument based on equilibrium liquids 
Q, however, suggests that the distribution should fall off 
differently for different potentials: the large force tail de- 
pends on the pair-distribution function at small r, which 
in equilibrium varies approximately as exp(— I/(r)/fcT). 
By this reasoning, one would predict a Gaussian tail for 
P{F) for a harmonic potential. Indeed, we do find a 
Gaussian tail except at T = near 00 (the peak of 
Pj{(j)c)), where both jammed and unjammed configura- 
tions can exist. In this region, we obtain a different re- 
sult for P{F), depending on whether we scale F by its 
average value before or after taking the configurational 
average. If, for each configuration, the forces F are scaled 
by the average value for that configuration, (F) , and the 
resulting distribution P{F/ {F)) is then averaged over all 
configurations, the result is as shown in Fig. ||(a). On the 
other hand, if the forces are scaled by the average over 
all configurations, {{F)), the resulting P{F/ {{F))) is as 
shown in Fig. ||(b). We see that a,t cj) = 0.644, near the 
peak of Pj{4>c) for N = 1024, the high force tail falls off 
much more slowly in this latter case than in the former, 
where each configuration is averaged separately. As the 
packing fraction is increased above (/)c, there is less differ- 
ence between the curves P{F/{{F))) and P{F/{F)). This 
trend can be explained since near (pc (F) varies dramat- 
ically from configuration to configuration. The relative 
size of fluctuations in (F) decreases with increasing (j), 
and the exponential decay of P{F) thus disappears. As 
the size of the system increases, the width of the region 
in (f) over which the high-force tail is exponential will de- 
crease. However, one can always tune (p close enough to 
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(j)o to observe the exponential tail. 

The shape of the tail of P{F/ {{F))) can be computed 
analytically, given a few simple assumptions. We have 
found that the average force in a configuration, (F), 
is proportional to the pressure p at T — near the 
onset of jamming. So from Fig. |l|(a), we know that 
(F) = Fo{(j) — (j)c) for the harmonic potential. We assume 
that the jamming threshold, 0c, is distributed as a Gaus- 
sian centered at (f>Q with width w as found in Fig. |^(b). 
Finally, we assume (as shown in Fig. ^(a)) that the tail 
of P{F/ (F)) for individual configurations is given by the 
equilibrium argument, which implies a Gaussian tail cen- 
tered at F/{F) = with width ap- (Here, Fq, ap, 0o, 
and w are parameters that can be obtained from the sim- 
ulation data.) In the large F limit, we find 

P(F)oc /^d</)c74re"'''^^'^''^'"^^e-(^=-*'')'/(2-') (1) 
Jo \P) 
^ eM-F/{{F))) 
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We have also studied Hertzian potentials and find re- 
sults similar to those for the harmonic potential shown 
in Fig. ^ and Eq. |^. In experimental granular sys- 
tems, where the interparticle potential is expected to be 
Hertzian at contact, P{F) has an exponential tail even 
for a single configuration We speculate that this is 
due to friction in the laboratory system, which allows het- 
erogeneities from region to region within a single sample. 

We have shown that for a finite-size system the jam- 
ming phase diagram looks somewhat different from the 
one sketched in the inset to Fig. |^(a). Instead of a well- 
defined point J, we find that there is a region of cen- 
tered around 0o with width w, in which both jammed and 
unjammed states can exist. As the size of the system in- 
creases, this region shrinks to the point J. Effects not 
included in our simulations, such as the presence of fric- 
tion, non-spherically symmetric potentials, or anisotropic 
packing (such as sequential packing under gravity) may 
prevent this region from disappearing. 

In some ways, point J in the phase diagram resem- 
bles a critical point: there is power-law scaling (Fig. |^), 
P{F/ {{F))) has a robust exponential tail independent 
of potential, and configurations are not self-averaging. 
In the context of foam, it has also been speculated that 
point J corresponds to rigidity percolation |p^ . However, 
near J the behavior differs from ordinary critical behav- 
ior, where configurations are not self-averaging once the 
correlation length exceeds the system size. There are no 
fiuctuations near J; that is, an unjammed (jammed) con- 
figuration will be unjammed (jammed) everywhere. This, 
as well as the fact that no bonds exist at packing fractions 
below J, makes this transition also different from rigidity 
percolation. Moreover, even though we find power-laws 
near J, there is a jump from Z = to Z = Zc and the 



exponents depend on the potential but not on dimension. 
Thus, point J has rather special properties. 

We thank J. -P. Bouchaud, P. Chaikin, D. Durian, G. 
Grest, D. Levine, J. Socolar, and T. Witten for helpful 
discussions. We acknowledge support from NSF-DMR- 
0089081 (CSO,SRN) and DMR-0087349 (CSO,AJL) and 
computing resources from the High Performance Com- 
puting Research Facility at Argonne National Lab. 



[1] See Jamming and Rheology ed. A. J. Liu and S. R. Nagel 
(Taylor & Francis, N. Y., 2001), and references therein. 

[2] A. J. Liu and S. R. Nagel, Nature 396 N6706, 21 (1998). 

[3] C. S. O'Hern, S. A. Langer, A. J. Liu, S. R. Nagel, Phys. 
Rev. Lett. 86, 111 (2001). 

[4] The peak in P{F) is more generally indicative of a yield 
stress, and appears also in crystals and ID chains. 

[5] D. M. Meuth, H. M. Jaeger, and S. R. Nagel, Phys. Rev. 
E 57, 3164 (1998); D. L. Blair, N. W. Meuggenburg, A. 
H. Marshall, H. M. Jaeger, and S. R. Nagel, Phys. Rev. 
E 63, 041304 (2001). 

[6] H. A. Makse, D. L. Johnson and L. M. Schwartz, Phys. 
Rev. Lett. 84, 4160 (2000). 

[7] F. Radjai, M. Jean, J. -J. Moreau, S. Roux, Phys. Rev. 
Lett. 77, 274 (1996); C. Thornton, KONA Powder and 
Particle 15, 103 (1997); M. L. Nguyen and S. N. Copper- 
smith, Phys. Rev. E 62, 5248 (2000). 

[8] F. H. Stillinger and T. A. Weber, Science 225, 983 
(1984). 

[9] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. 

Vetterling, Numerical Recipes in Fortran 77 (Cambridge 

University Press, New York, 1986). 
[10] R. J. Speedy, J. Chem. Phys. 110, 4559 (1999). 
[11] D. N. Perera and P. Harrowell, Phys. Rev. E 59, 5721 

(1999). 

[12] D. J. Durian, Phys. Rev. Lett. 75, 4780 (1995); Phys. 
Rev. E 55, 1739 (1997). 

[13] T. G. Mason, M.-D. Lacasse, G. S. Grest, D. Levine, J. 
Bibette, and D. A. Weitz, Phys. Rev. E 56, 3150 (1997). 

[14] This result suggests that states at <j)c are not "fragile", 
as in M. E. Gates, J. P. Wittmer, J.-P. Bouchaud, and P. 
Claudin, Phys. Rev. Lett. 81, 1841 (1998). It is possible 
that states below (j>c can jam when a shear stress is ap- 
plied. Such states can support a pressure but are fragile 
because they will yield when the shear stress is reversed. 

[15] J. G. Berryman, Phys. Rev. A 27, 1053 (1983); A. S. 
Clarke and J. D. Wiley, Phys. Rev. B 35, 7350 (1987). 

[16] S. Torquato, T. M. Truskett, and P. G. Debenedetti, 
Phys. Rev. Lett. 84, 2064 (2000). 

[17] J. D. Bernal, Nature 188, 910 (1961); J. D. Bernal, J. 
Mason, and K. R. Knight, Nature 194, 956 (1962). 

[18] G. Y. Onoda and E. G. Liniger, Phys. Rev. Lett. 64, 
2727 (1990). 

[19] F. Bohon and D. Weaire, Phys. Rev. Lett. 65, 3449 
(1990); D. J. Jacobs and M. F. Thorpe, Phys. Rev. E 
53, 3682 (1996). 



4 



